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Abstract. - We study the response of a two-dimensional hexagonal packing of rigid, friction- 
less spherical grains due to a vertically downward point force on a single grain at the top layer. 
We use a statistical approach, where each configuration of the contact forces is equally likely. 
We show that this problem is equivalent to a correlated g-model. We find that the response 
displays two peaks which lie precisely on the downward lattice directions emanating from the 
point of application of the force. With increasing depth, the magnitude of the peaks decreases, 
and a central peak develops. On the bottom of the pile, only the middle peak persists. The 
response of different system sizes exhibits self-similarity. 



Force transmissions in (static) granular packings have attracted a lot of attention in recent 
years [1-13]. Granular packings are assemblies of macroscopic particles that interact only via 
mechanical repulsion effected through physical contacts. Experimental and numerical studies 
of these systems have identified two main characteristics. First, large fluctuations are found 
to occur in the magnitudes of inter-grain forces, implying that the probability distribution 
of the force magnitudes is rather broad [4]. Secondly, the average propagation of forces - 
studied via the response to a single external force — is strongly dependent on the underlying 
contact geometry [3,11-13]. 

The available theoretical models capture either one or the other of these two aspects. 
The scalar q-model [5] reproduces reasonably well the observed force distribution, but yields 
diffusive propagation of forces, in conflict with experiments [11, 12]. Continuum elastic and 
elasto-plastic theories [6] predict responses in qualitative agreement with experiments [7-10], 
but they provide a description only at the average macroscopic level. More ad-hoc "stress- 
only" models [2] include structural randomness, but its consequences on the distribution of 
forces are unclear. In other words, an approach that produces both realistic fluctuations and 
propagation of forces in granular materials from the same set of fundamental principles is still 
called for. 

A simple conjecture, which could provide such a fundamental principle for all problems of 
granular statics, has been put forward by Edwards years ago [14, 15]. The idea is to consider 
all "jammed" configurations equally probable. A priori, there is no justification for such an 
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ergodic hypothesis, but its application to models of jamming and compaction has been rather 
successful [16]. Its extension to the forces in granular packings is in principle straightfor- 
ward: sets of forces belonging to all mechanically stable configurations have equal probability. 
However, in an ensemble of stable granular packings, two levels of randomness are generally 
present [2]. First, the force geometry clearly depends on the underlying geometrical contact 
network, which is different in different packings. Secondly, randomness in the values of the 
forces is present even in a fixed contact network, since the forces are not necessarily uniquely 
determined from the contact network. Instead of considering both levels of randomness si- 
multaneously, a natural first step is thus to obtain the averages for a fixed contact geometry, 
and then possibly to average over the contact geometries. 




Fig. 1 - (a) The model: (N + 2) xp array of hexagonally close-packed rigid frictionless spherical grains 
in two-dimensions (drawn for odd N). At the top, there is only a single vertically downward point 
force applied on one particle. At the boundaries, little gray circles appear on interfaces where the 
contact forces are non-zero, (b-c) Schematically shown forces on the jth grain in the ith layer: (b) 
i < N, (c) i = (N + 1), the bottom reaction Wn+i is shifted upwards for clarity; Fm^ > Vm. 



While such a method has recently been shown to produce single inter-grain force probability 
distributions in fixed geometry that compare well with experiments [17], in this Letter we 
demonstrate that it also leads to an average response function qualitatively in agreement with 
experiments. More precisely, we determine the behaviour of the response of a two-dimensional 
hexagonal packing of rigid, frictionless spherical grains placed between two vertical walls (see 
FigDi due to a vertically downward force F applied on a single grain at the top layer. 
Experimentally [11, 12] it was found that a force F applied to the top of a hexagonal packing 
of photo-elastic particles propagates mainly along the two downward lattice directions. We 
define the response of the packing as (Wij) — (W^j ) /F where Wij and are the vertical 

force transmitted by the (i,j)th grain to the layer below it respectively with and without 
the external force F, and the angular brackets denote averaging over all configurations of 
mechanically stable contact forces with equal probability. 

To start with, we describe a method for assigning the uniform probability measure on the 
ensemble £ of stable repulsive contact forces pertaining to a fixed geometrical configuration 
of P rigid, frictionless two-dimensional grains of arbitrary shapes and sizes (for a rigorous 
geometrical description of a granular packing, see Ref. [15]). The directions of the forces 
are fixed at each of the Q contact points, and one can represent any force configuration 
by a column vector F consisting of Q non-negative scalars {Pfc} (with k = 1, . . . , Q) as its 
individual elements. These elements satisfy 3P Newton's equations (3 equations per grain: 
two for balancing forces in the x and y directions and one for balancing torque), which can 
be represented as A • F = F ex t- Here, A is a 3P x Q matrix, and V ext is a 3P-dimensional 
column vector representing the external forces. If we assume 3P < Q [20], then there is no 
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unique solution for F. Instead, there exists a whole set of solutions that can be constructed via 
the three following steps: (1) one first identifies an orthonormal basis {F^"} (1 = 1,..., cLk = 
Q — 3P) that spans the space of Ker(A); (2) one then determines a unique solution F^ ^ of 
A • F^ ) = F ex t by requiring F^.F^ = for I = 1, . . . ,c?k; and (3) one finally obtains all 

Q-3P 

solutions of A ■ F = F ex t as F = F' ) + ft ^ > where /z are real numbers. This implies 

i=i 

that £ is parametrized by the /; 's belonging to a set S obeying the non- negativity conditions 
for all forces. The uniform measure on £, which is usually compact [23], is thus equivalent to 
the uniform measure d[i = Y[ dF^ S(A ■ F — F ext ) <d(Fk) = Y[ dfi on S. 

k I 

In our model, the grains are spherical, so that the dimension of A reduces to 2P x Q. We 
consider the force F and the weights of the individual grains as the non-zero elements of F ex t , 
while F is composed of all inter-particle and non-zero boundary forces. Simple counting then 
shows that Q — 3Np + 5p + N + 2 (see Fig. QJ. The matrix A represents two equations per 
particle (see Fig. [Tj (b-c)) [20] 
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i.e., 2(N + 2)p equations all together, implying that dx = N + 2+(N + l)p. We choose F^'^'s 
for i < (N+l) and F^'^'s for i < N, 1 < j < N to parameterize £ . Once these forces are fixed, 
all the others are uniquely determined by solving Eq. Q layer by layer from top down [24]. 
It is easily seen that the number of these parameters is indeed <1k, as it should be. Clearly, 
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implying that the set 5' of allowed values of Gij's is compact. However, since the non- 

(i i) 

negativity conditions for F 3 
model is actually unbounded. 



and F^'^'s provide only lower bounds for F^' 1 '^, S in this 



The remedy we use is to fix the F3 values: indeed, as can be seen in Eq. QJ, the values 
of the Wij depend only the Gjj's so that in this model the precise values of if 1 -* have no 
physical meaning. Nevertheless, one has to be careful: notice that the Gjj 's are differences 
of the physical contact forces and thus they are allowed to become negative in magnitude. In 
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Fig. 2 - Colour plots for N = 35 and m — 0: (a) mean response; (b) the standard deviation of the 
response. 



fact, Eqs. and Q together imply that if F^' 1 ' < 2[F + (i — l)mg]/y/3, then the positivity 

(i i) 

requirement of the _F 4 might further restrict the choice of Gjj values within the bounds 
of inequality (J3J. In this Letter, we fix the mag nitudes F 3 (4,1) = 2[F + (i - l)mg]/V3 = F 
so that all values of dj within the bounds of inequality © are allowed (details of the cases 
^ 3 (i,1) < F appear elsewhere [18,19]) This arrangement reduces the uniform measure over £ 
to the uniform measure on S' , which is a (N + l)]3-dimensional polyhedron. 

To evaluate (Wij) = -r? [ Wij TTdG^j, where N = J s , FJ < . dGij is the normalization 

J S' i„i 

constant, we define 



q itj = \V3(Gij + F^ ) )/2 + mg/2] /W itj , 



(4) 



where qij is the fraction of Wij that the (i, j)th particle transmits to the layer below it toward 



the left, i.e., Ff 3 ' = 2q t jW t j / V?> and Ff 3 ' = 2(1 - q l j)W l ^ /V3. Equation <gj) then reduces 
Eq. © to < qi.j < 1. Clearly, Woj are the external forces applied on the top layer. For 
i > 0, Wi j is a function of q^i for k < i, since 



Wij = (1 - Wi_ij_i + qi-x,j Wi-ij + mg. 



(5) 





-0.4 -0.2 0.2 0.4 

x/z 

Fig. 3 - Behavior of (W(i,j)) for m = in reduced co-ordinates x and z: (a) scaling of (W(x,z)} 
with system size for \x\ < z/2 at two z values (for clarity, z = 0.8 and z = 0.6 have been shifted 
upwards by one and two units respectively); (b) data collapse for (W(x, z)) \ x=z /2 at three N values. 
See text for further details. 
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It may seem from Eq. ^ that in the hexagonal geometry of Fig. ^ one simply recovers 
the g-model [5]. There is however an important subtlety to take notice of. In the g-model, the 
g's corresponding to different grains are usually uncorrelated, while in our case, the uniform 
measure on S' implies, from Eq. Q), that 

I] dGij = 2( w +Dp Yldqij Wij (q)/^ N+1 ^ 2 . (6) 

id i,3 

Due to the presence of the Jacobian on the right hand side of Eq. @) , the uniform measure 
on S' translates to a non- uniform measure on the [N + l)p-dimensional unit cube formed by 
the accessible values of the g's. 

Notice an important artifact of this approach: the joint probability distribution Y[ Wij(q) 

depends on the qij values over the whole system, thereby making cfoj's correlated with each 
other. In fact, the induced probability P(qij) for a single grain does not only depend on the 
number of layers present above the grain, but also on the number of layers below it. It is 
thus clear that the forces in this model do not propagate top down, as they do in hyperbolic 
"stress-only" models [2]. 

For massless grains (m = 0), it is clear that (i) for F — 0, W^ = 0V(i, j) (ii) the (Wij) 
values scale linearly with F (hence, we use F — 1), and (hi) (Wij) = outside the triangle 
formed by the two downward lattice directions emanating from = (0,jo)j the point of 
application of F. The (Wi j)'s, evaluated numerically via the Metropolis algorithm on these 
g's, appear in Fig. 

Our simulation results for (Wi.j) and the standard deviation SWij — J (Wf^) — (Wij) 2 
within the triangle are plotted in Fig. [21 using the built-in cubic interpolation function of 
Mathematica. Outside the triangle, (Wij) = appears in deep indigo; the largest (Wij) value 
within the triangle appears in dark red; and any other non-zero (Wij) value is represented 
by a linear wavelength scale in between [25]. We find ViV that (a) the (Wij) values display 
two single- grain- diameter-wide symmetric peaks that lie precisely on the two downward lattice 
directions emanating from jo, (b) the magnitudes of these peaks decrease with depth, and (c) 
only a central maximum for (Wij) is seen at the very bottom layer (i = N+l). The standard 
deviation 5Wij has a similar shape to (Wi j), although the peaks are less pronounced. 

We further define x — (j — jo)/(N + 1) and (j — jo + 1/2) /(N + 1) respectively for 
even and odd i, and z = i/(N + 1) [see Fig. ^ in order to put the vertices of the triangle 
formed by the locations of non-zero (Wij) values on (0,0), (—1/2, 1) and (1/2, 1) ViV. The 
excellent data collapse shown in Fig. [3] indicates that the (Wij) values for |x| < z/2 scale 
with the inverse system size [Fig. Efa); we however show only three z values], while the 
(Wij) values for ~ z/2 lie on the same curve for all system sizes [Fig. EIb)]. The data 
suggest that in the thermodynamic limit N — > oo, the response field (W(x, z)) scales ~ \/N 
for |x| < z/2, but reaches a non-zero limiting value on |a;| — z\/z < 1. We thus expect 
lim (W(x, z)) | | a .|_ z ^ 2 > (W(x, z)) 1 1 2 ,| <2 ^ 2 ^ z < li or equivalently, a double-peaked response 

field at all depths z < 1 for large N. 

We have not found a simple explanation for such scaling behaviour of (W(x, z)). It however 
turns out that exact analytical expressions can be obtained for all moments of Wij V(i,j), 
for any N. The detailed calculations appear in Ref. [18]. 

In view of the self-similarity of (W(x, z)) that we observe for different system sizes in Figs. 
OJb-c), it seems natural that we also study the same properties for m ^ 0. In this case, 
(W$) + and (Wij) 9^ F. To minimize the effect of the boundaries in the regions around 
j = jo, w e have used p — 2N + 5. For m ^ 0, the relevant scale for the magnitude of F is 
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Fig. 4 - Scaling properties of (w(x,z)), analogous to Figs, fflb-c). for /3 = 100 and three A^ values 
(a.u. = arbitrary units). 



obviously a — F/mg. For q ^ 0, just like in the case of m — 0, we observe a double-peaked 
response, and the peaks are still single grain diameter wide. Furthermore, for a given value of 
N and increasing a, the magnitude of the response on x — z/2 decays more slowly, i.e. the 
peaks penetrate the packing to progressively higher values of z. In order to avoid repetition, 
we do not use colour figures like Fig. EJa) to demonstrate this behaviour, but the trend of the 
data clearly indicates that for fixed N, one should recover the results corresponding to m = 
in the limit a — > oo. 

It is clear from the qualitative behaviour described in the above paragraph that in order 
to obtain scaling with increasing N, one also needs to scale a in some way. To this end, we 
define (3 — a/(N + 1) and keep j3 constant for increasing N. The corresponding graphs are 
shown in Fig. 0]for j3 = 100. The fact that the self-similarity in Fig. 0]for different system 
sizes is not as striking as in Figs. EJb-c) suggests that there is more to the story of scaling 
properties. It is likely that the full scaling properties can be unraveled only at much higher 
values of N, but unfortunately, simulations with TV values significantly higher than 50 require 
impractically long times. 

In summary, we find that assigning equal probability to all mechanically stable force con- 
figurations for rigid, frictionless spherical grains (with or without mass) in a two-dimensional 
hexagonally close-packed geometry yields a double-peaked response. The peaks are single 
grain diameter wide, they lie on the two downward lattice directions emanating from the 
point of application of F. With increasing depth, the magnitude of the peaks decreases, and 
a third peak starts to develop directly below the applied force. Near (and on) the bottom 
layer only the middle peak persists; i.e., the response becomes single-peaked. As the number 
of layers is increased, the transition from double to single peak takes place deeper in the 
packing. Moreover, for grains each with mass m, the peaks penetrate the packing deeper with 
larger F . The standard deviation of the response is similar in shape to the response, but the 
peaks are weaker. 

We emphasize that the results presented here are obtained for the boundary condition 
^ 3 (i,1) > F . The case Fg*' 1 ^ < Fq and other kinds of boundary conditions have been analyzed 
elsewhere [18, 19]. These results together indicate that the quantitative behaviour of the 
response depends crucially on the side forces (i.e. boundary conditions) — this feature is 
consistent with other theoretical approaches [7]. In particular, we note that the transition to 
a single-peaked response does not take place for F^ 1 ' 1 ^ sufficiently small. 

We also note that the double-peaked structure of both the mean response and the stan- 
dard deviation of the response is in qualitative agreement with experiments [11-13], but the 
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fluctuations observed in the model are much weaker than found in experiments [12]. Another 
crucial difference between our and the experimental results is that in this model the peaks are 
single-diameter wide, while in experiments the peaks widen with depth [12]. This difference 
probably stems from the fact that in experiments the effect of inter-grain friction can never be 
neglected. Presence of friction would also certainly make fluctuations in the response SWij 
stronger. A study of the effects of friction on the response along the lines of [26] is therefore 
an important direction for future work. 

It is a pleasure to thank J. -P. Bouchaud, D. Dhar, J. M. J. van Leeuwen, B. Nienhuis, 
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